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SUMMARY 

Using measured values of two-dimensional forces in a magnetic actuator, equations of 
motion for an active magnetic bearing are presented. The presence of geometric coupling between 
coordinate directions causes the equations of motion to be nonlinear. Two methods are used to 
examine the unbalance response of the system: simulation by direct integration in time: and 
determination of approximate steady state solutions by harmonic balance. For relatively large 
values ot the derivative control coefficient, the system behaves in an essentially linear manner, but 
for lower values of this parameter, or for higher values of the coupling coefficient, the response 
shows a split ot amplitudes in the two principal directions. This bifurcation is sensitive to initial 
conditions. The harmonic balance solution shows that the separation of amplitudes actually 
corresponds to a change in stability of multiple coexisting solutions. 

INTRODUCTION 

In the short time since the introduction of practical magnetic levitation for rotating shafts, 
significant progress has been made in designing and modelling the performance of magnetic 
bearings, to the point where they are practical for a variety of applications. These include machine 
tool spindles, pumps, compressors, gyroscopes and momentum wheels [1 ]. Some applications, 
those with signiticant rotor flexibility and coordinate cross coupling, present difficult challenges, 
however. In this category fall some pumps and many compressors and turbines, including gas 
turbines. Recent papers, such as those of Williams, et al. [2], Lee and Kim [3], and Nonami , et 
al.[4] have confronted limitations of vibration control at high speeds for flexible rotors using 
magnetic bearings. Rotor dynamic stability and robust vibration control in these machines will 
increasingly depend on an understanding of the system dynamics that goes beyond traditional 
linear models. 


This paper presents equations of motion and simulations of a two-degree-of-freedom 
system subject to forces from an actively controlled magnetic actuator, where the control is linear 
but the forces from the actuator include coordinate coupling of a form found in experimental 
measurements. The resulting equations are nonlinear and coordinate -coupled, and the response 
contains several features found only in nonlinear systems. 

ACTUATOR FORCES 

A schematic of an active magnetic bearing is shown in Figure I. It consists of two 
opposed pairs ot electromagnets arranged around a shaft. Each of the magnet pairs will be 
controlled independently. The control is linear, taking as input the shaft displacement along the 
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axis of the magnet pair, and producing as output a variation ot the magnetic flux B. This type ot 
control is an idealization from the practical case, where control is more easily exerted on either the 
current or the voltage applied to the magnet coils. 



Figure 1 . Magnetic bearing actuator. Sensors and control not shown. 

The force between a magnet and a ferromagnetic part is the negative ot the derivative ot 
magnetic field energy with respect to virtual displacement ot the pail. The traditional model toi 
forces from a magnetic actuator in a magnetic bearing is based on one-dimensional magnetic circuit 
theory in which it is assumed that the lines ot magnetic flux cross the air gaps ot the beanng in 
straight lines It is generally assumed that the direction of the force is along the axis ot symmetry 
of the magnet. In a previous paper [5] measurements were described in which the forces from a 
ma°net with curved pole faces acting on a circular shatt had both an attractive component and a 
component normal to the axis of the magnet when the shaft was given a displacement away trom 
this symmetry axis. The ratio of forces was found approximately proportional to the normal 
coordinate. 

In the calculations that follow, both principal forces and normal forces are assumed to act 
The principal force from each magnet is modelled by this one-dimensional circuit theory, in which 


where a = pole cross section area 

B = magnetic flux 

|i 0 = permeability of free space 

Often the flux is written as a function of the coil current i. wire turns N and gap h. leading to 

c MolNifa , 
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In the present work, the controlled quantity is assumed to be the flux B, so Equation ( I ) 
will be used to model the attractive force between each magnet and the shaft 

Based on the measurements of [5], the normal force is considered to be 

F n = ocF p (2) 

where the coupling parameter a is empirical. The directions p,n correspond to directions x,y or 
y,x as appropriate for each magnet. 


Linear Flux Control 

Consider the case in which a point mass is acted on by two opposed magnet pairs, shown 
schematically in Figure I . In an attempt to linearize the system so that the mass effectively is 
subject to a linear restoring force like a spring stiffness, a large bias flux B b is introduced in each 
magnet gap, and a control flux B c is superposed. If the bias flux is equal in both magnets of the 
vertically oriented pair 


B| b = B 3b = B b (3) 

This constitutes the special case in which the bearing is not required to support a steady load, as in 
a vertically oriented machine. 

It the control flux is made equal in magnitude but opposite in sign in the two magnets 

B| C = -B 3c = B c (4) 


then the resultant force in the y direction 

By (1.3) = [(B i b + B i c )~ - (B 3b + B 3c )-] 


(5) 


is linear in the control flux 


F y (1.3) 


4a 


B b B c 


( 6 ) 


Coordinate Coupling 

The approach above ignores the normal force found by experiment and calculation. This 
normal force results in the need for control in the x direction as well. Suppose that an additional 
magnet pair is placed on the x axis, subject to the same bias and control flux as the y-direction pair. 
For instance, let the control flux in each magnet be proportional only to the distance from that 
magnet, neglecting for the present any derivative control. That is 

B c = - KB b x p (7) 

where the subscript p denotes the principal axis, the axis of symmetry of the magnet. 


355 



Then, adding all the forces in the x direction 


Fx = i5lEf K x-e.x(i + KV)l < 8 > 

Uo L 2 J 

Similarly, the y direction force is 

F y = ^~b|Ky-^y(l +K 2 x 2 )] (9) 

It is seen that the coordinate coupling embodied in the parameter cc has two effects, it 
obviously yields a force normal to the displacement, but it also attenuates the lestoiing foice, 
reducing the effective stiffness. Both these effects occur even though the independent-axis control 
algorithm is assumed perfect in its ability to give a linear stiffness chaiacteristic in the piincipal 
direction. 

The forces given by Equations (8) and (9) are derivable from a potential energy function 

v = 4aBgr/jc _a\( x 2 + y 2) . Q, K 2 x 2 y 2 l < 10 > 

Ho L'2 4 ' ' J 4 j 


This quantity may be rendered dimensionless 

v = (K . A.) (x 2 + Y 2 ) - ^ K 2 X 2 Y 2 < 10a) 

using the parameters 

K = kc A = ac 

X = x/c Y = y/c 

(ID 

V = v — 

4aBj;c 

Figure 2 contains plots of the dimensionless potential energy at different levels of the 

coupling parameter A [6], The attenuation of principal stiffness with increasing A is evidenced by 
the decrease in amplitude of the potential, and the normal forces are evidenced by the tendency of 
the potential surface to sag along lines at 45° to the principal axes. One way to considei the 
implications of these plots is to visualize them as solid surfaces upon which a mass in the shape of 
a small sphere might roll if given an initial velocity and/or displacement from the origin. This is a 
physical analog to the free vibration of a mass subject to the force equations given above. 
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FORCED RESPONSE 


If a derivative feedback is added to the proportional feedback 

Be = ■ K-BfaXp - 


the force equations 

F _ 2aBg r 2Kx + ax (j + K 2 y 2)] _ 2yx + 2aKyxyy + crrxy 2 (13) 

Ho 

2 

p 2 Ky + ay(l + K 2 x 2 )] - 2yy + 2aKyyxx + airyx 2 (14) 

y Ho 

show a coupling of the proportional and derivative terms in each coordinate, as well as a coupling 
of both these terms to the other position coordinate. 

Consider a mass m subject to the above actuator forces in addition to an unbalance fore in 
function. Newton’s second law is written for each of the coordinate directions, with the resultin 
equation of motion in x as 


m x = d5®b[ K x-^x(l + K 2 y 2 ) + YX - aKyxyy - ktrxy 2 

Ho L 2 


+ meocrcosoot (15) 


or 


rnx = ~ 4at3 i [kx + 7X - £x ( 1 + K 2 y 2 + otKyyy + W^ 2 )] + meio 2 cos(i>t ( 1 6) 
Ho L 21 “ 


and in y 


my = 


-4aBi: 


Ho 


Ky + yy - % 1 1 + k 2 x 2 + aKyxx + 


IT 2 


+ meorsintot 


17) 


By choosing the additional nondimensional parameters 

K = k/c E = e/c 


T = 


t(D n 


r = Ycco n 


Q = co/to n 


®n = 


1 4aKB^ 

V m|i 0 


(18) 


the system of equations can be nondimensionalized. Here, (On is the natural frequency the system 
would have if there were no coupling present (a = 0). As Equation (10) indicates, for any positive 
value of a, the actual natural frequency will be reduced from this value. 
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The nondimensional forms of the equations are 

x = " j^"( xx + EX' ' ^ x ( 1 +K“Y~ + 2KrYY' + r~Y ' 2 )j + ELTcosQT (19) 

Y " = ' K ( KY + rY ' ' 1 +K2X2 + 2Krxx ' + r2x ' 2 )) + ECTsinQT (20) 

where 1 and 11 indicate differentiation with respect to dimensionless time T. These equations may be 
integrated in time after assigning values to the system parameters K, T, A and E. along with 
appropriate initial conditions for all of the quantities X, Y. X’ and Y'. It should also be noted that 
the toms of the forcing function in the final terms of Equations ( 19) and (20) also constitute initial 
conditions, in the form of an assumed phase angle for the forcing function at T = 0. 

ANALYSIS 

Two methods are used to examine the response of the system to unbalance forcing: 
numerical integration from arbitrary initial conditions using a fourth-order Runge-Kutta algorithm: 
and approximate calculation of steady state solutions by the harmonic balance method. 

Natural Frequency 

Because the potential energy of the system as shown in Figure 2 is not described by a 
surface ot revolution, the natural frequency of the system depends on the initial conditions. If 
initial conditions of finite displacements and zero velocities are chosen, the period of oscillation in 
tree vibration is a function ot the displacements, or alternatively, of the radius and angle of the 
initial conditions. Figure 3 indicates this dependence. The results were obtained by direct 
numerical integration, and dissipation was not included. Part (a) shows the period as a function of 
radius for a fixed angle from the x-axis (22.5°) and part (b) shows the effect of the angle of the 
initial condition for a fixed radius (0.5). The two periods both increase with radius, a characteristic 
of systems with softening stiffness. The period in Y is different from that in X because the path of 
the oscillation does not in general pass through the origin. There is a single period only when the 
initial conditions lie along a line at 45° to the two axes, as shown in Figure 3b. 




(a) Fixed angle from x axis (b) Fixed radius 

Figure 3. Effect of Initial Conditions on Natural Period 
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Forced Response by Numerical Integration 


When the derivative control coefficient T is made finite so the system dissipates energy, the 
response of the system to unbalance forcing can be calculated. The ettects of thiee paiameteis aie 

examined: K, T. and A. 

The proportional control coefficient K determines how rapidly the flux, and thereby the 
force, from a magnet is reduced as the shaft approaches that magnet. A value ot K=1 would cause 
the flux to be reduced to zero at contact (not counting the contribution from the derivative control 

coefficient T). Larger values of K would correspond to "stiffer" bearings. Because of the form ot 

nondimensionalization ot Equations ( 19) and (20). however, an increase in K while holding T 
constant causes a decrease in the effective dissipation coetlicient. by virtue ot a change in the 
natural frequency. This must be considered when interpreting the results of parametric studies. 

since a straightforward increase in the dimensional quantity k ( 1/m) would not attect the 
dissipation, or derivative control, coefficient. The tact that K cannot be eliminated from the 
equations of motion is a result of the essential nature of nonlinear systems. 

The measurements of [5] indicate that A=0. 1 5 is a reasonable value tor the cooidinate 
coupling coefficient A. In the calculations below, A is varied from 0.05 to 0.25.The values ot T 
were chosen to provide dissipation of the same order as in a linear system having damping ratios 
between 0.1 and 0.3 . In all the results presented here, the unbalance magnitude is E=0.l . 

Studies examining large ranges ot parameter combinations are planned but are beyond 
the scope of the present work. Nevertheless, much can be learned from a limited parametric 
study. 

Figure 4 shows the effect of increasing the coefficient K from 1 .0 in part a to 3.0 in part b 
to 5.0 in part c. As noted above, increasing K alone results in a smaller value ot the dei ivative 
control coefficient. With this in mind. Figure 4 still indicates an important feature of the system: 
that for some combinations of parameters, the response exhibits a split, with the motion on one 
axis having a much higher amplitude than that on the other axis. Associated with this split is a 
sudden jump in one of the amplitudes as the frequency is increased. In fact, one of the solutions ot 
Figure 4c extends beyond an eccentricity of 1 .0, which in the physical case would result in solid 
contact. At some frequencies near the natural frequency, however, two solutions exist that are 
both within the physical bounds of the system. Furthermore, the solutions are dependent on the 
initial conditions. For the case shown, the integration begins with both shaft position and velocity 
equal to zero. The numerical integration proceeds until all transients have decayed and the peak 
amplitudes in the two directions are sampled. The forcing function, the final terms in Equations 
(19) and (20), also imposes an implicit initial condition by virtue of its assumed phase. In tact, the 
cosine portion of the forcing function begins with a step imposition ot force at time t=0. although 
all transients associated with this discontinuity have decayed before the amplitudes are sampled. 

If. however, the cosine and sine parts of the force are exchanged, the solutions toi X and Y aie 
also found to have exchanged places. This is in marked contrast to a lineai system, wheie the 
amplitude is unique after transients have decayed. 
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Figure 4. Effect of variation in dimensionless proportional control coefficient K. 




The effect of the coupling parameter can also be examined. Figure 5 illustrates the effect ot 
increasing the value of A, while K is held constant at 0.3. The values ot all parameteis except A 
are equal to those of the case shown in Figure 4b. As A is increased beyond a threshold value, 
between 0.15 and 0.25. multiple solutions appear near the natural frequency. In this set ot plots, 
the natural frequency is a constant, making this parametric variation somewhat easier to interpret 

than the previous one. 

Reduction of the derivative control coefficient T can also bring about a situation with 
multiple solutions, as shown in Figure 6, as can an increase of the unbalance eccentricity, not 
shown. 

Thus, the bifurcation seems mostly to be an amplitude-driven phenomenon, such that when 
a critical amplitude is exceeded, the solutions split. In all cases, the split is initial-condition- 
dependent. 

Solution by Harmonic Balance 

Another approach to examining the steady-state response of a nonlinear system is the 
harmonic balance method, which is approximate but analytical rather than numerical. It has the 
advantage that both stable and unstable solutions can be located, whereas numerical integration can 
locate only stable solutions. 

The method consists of assuming steady solutions of the form 


X = C cosQT + DsinQT 
Y = G cosQT + HsinQT 


( 21 ) 


( 22 ) 


where C D G, and H are to be determined. Equations (2 1 ) and (22) are differentiated and 
substituted into the equations of motion. The resulting powers of trigonometric functions are 
expanded using trig identities, after which the harmonics higher than 1 are neglected. Because the 
truncation of higher harmonics is not performed until after the poweis ot trig functions aie 
expanded, the solution retains its nonlinear character, although the equations have been 
approximated. The resulting four algebraic equations for the constants C, D, G. and H are coupled 
and highly nonlinear and must themselves be solved by a numerical Newton-Raphson iteration LoJ. 
When die constants are found, the steady amplitudes can be calculated readily. 

Figure 7 shows the amplitudes obtained by harmonic balance for the case corresponding to 
Figure 5a° These results indicate that in the neighborhood of the natural frequency, four solutions 
actually exist. (Two are identical.) Two of the solutions are apparently unstable, but the harmonic 
balance method does not yield stability characteristics. Based on the results of numerical 
integration, however, it appears that the solutions corresponding to equal amplitudes for X and Y 
are unstable when they lie between the unequal solutions. Thus the jump in one ot the amplitudes 
stems from a change in that solution's stability. Where the equal-amplitude solutions lie below the 
unequal solutions, they are the stable ones. The unequal solutions are believed to exist at all 
frequencies, but are difficult to locate by Newton-Raphson beyond the range that is shown. 


362 










Amplitude Amplitude 


1 1.5 

Frequency ratio, ay(O n 


Frequency ratio, oV(D n 


1 1.5 

Frequency ratio, oV(j) n 


K = 3.0 
A = 0.15 

T = 0.2 



K = 3.0 
A = 0.15 

T = 0.6 


Figure 6. Effect of variation in derivative control coefficient f 







Figure 7. Multiple coexisting solutions by harmonic balance method. 

The close correspondence between the numerical and analytical results supports the validity 
ot both methods. Neither method alone is sufficient for a complete understanding, however, 
because the numerical solutions are dependent on initial conditions, and the analytical solutions 
provide no information on stability. 

The numerical integration can in fact be used to track the unstable solutions to a very limited 
extent by caretul choice of initial conditions, as shown in Figure 9 by the parts of the curves 
labelled "alternate solution". These initial conditions are based on the harmonic balance results, 
and tend further to support the validity of both methods. 



Figure 8. Limited unstable solution-following by numerical integration. 
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CONCLUSIONS 


Equations of motion and limited parametric studies are presented for the case of a magnetic 
bearing subject to flux control, with geometric coordinate coupling. There are two effects of the 
coupling parameter on the system potential energy: a reduction of the principal stiffness, and 
introduction of a nonlinear normal stiffness. 

The equations of motion are nonlinear and exhibit behavior that is distinctly different from 
that of linear systems. In forced response, the amplitudes of the system show bifurcations that are 
the result of changes in stability of multiple coexisting solutions. The stability seems mostly to be 

amplitude dependent, and the critical amplitude is a function of several parameters: K. T. and A. 

In the long term, successful implementation of magnetic bearings where large eccentricities 
may be encountered will depend on a deeper understanding of the nonlinear characteristics of the 
combined rotor-actuator-control system. 
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